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The 2nd-order upwind inviscid flux scheme implemented in the multi-block, 
structured grid, cell centered, finite volume, high-speed reacting flow code 
VULCAN has been modified to reduce numerical dissipation. This modification 
was motivated by the desire to improve the codes ability to perform large eddy 
simulations. The reduction in dissipation was accomplished through a 
hybridization of non-dissipative and dissipative discontinuity-capturing advection 
schemes that reduces numerical dissipation while maintaining the ability to capture 
shocks. A methodology for constructing hybrid-advection schemes that blends non- 
dissipative fluxes consisting of linear combinations of divergence and product rule 
forms discretized using 4th-order symmetric operators, with dissipative, 3rd- or 4th 
-order reconstruction based upwind flux schemes was developed and implemented. 
A series of benchmark problems with increasing spatial and fluid dynamical 
complexity were utilized to examine the ability of the candidate schemes to resolve 
and propagate structures typical of turbulent flow, their discontinuity capturing 
capability and their robustness. A realistic geometry typical of a high-speed 
propulsion system flowpath was computed using the most promising of the 
examined schemes and was compared with available experimental data to 
demonstrate simulation fidelity. 


Nomenclature 


F a = alternative split flux forms 

F c = product rule flux 

a = flux combination coefficient 

p = density 
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= Cartesian velocity components 
= static pressure 
= Static temperature 
= internal energy 
= Total energy 
= time 

= primitive variable 
= conserved variable 
= approximate flux jacobian 
= left state 
= right state 

= MUSCL scheme coefficient 
= pressure limiter coefficient 
= pressure limiter scaling coefficient 
= WENO normalized weight 
= WENO weight 
= WENO target weight 
= WENO smoothness indicators 
= WENO smoothness indicator limiter 
= discontinuity sensor coefficient 
= Ducros sensor constant to prevent division by zero 
= Larsson discontinuity sensor coefficient 1 
= Larsson discontinuity sensor coefficient 2 
= Larsson discontinuity sensor function 
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= cell length scale 
= ratio of the specific heats 
= eddy viscosity 

= van Driest wall damping coefficient 
= Smagorinsky constant 
= Smagorinsky grid filter width 
= rate of strain tensor 
= non-dimensional distance to the wall 
= wall shear velocity 
= dynamic viscosity 

= time interval normalized by the characteristic flow through time 
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I. Introduction 


N ASA research and development of air-breathing high-speed propulsion systems is currently being re-scoped. 

Consequently, investment in high-speed propulsion technology is being carefully examined for it's relevance 
and utility. However, research in computational fluid dynamics (CFD) tool development for hypersonic air-breathing 
propulsion has remained viable. This is due to the recognition that existing sub-scale demonstrator vehicles must be 
scaled up by a factor of 10 in order to satisfy real mission requirements. Unfortunately, such a 10X scaling cannot be 
easily accommodated by existing experimental facilities. Fortunately, sub-scale testing combined with CFD analysis at 
the 10X scale is considered a viable alternative if the ability of CFD to simulate the physical phenomena that 
contribute to off-design operability can be demonstrated. These phenomena are especially challenging in the case of 
the scramjet flowpath. 

The internal high-speed flows typically found in scramjet flowpaths pose a particular challenge for CFD 
because of complex, non-linear interactions between viscous wall bounded flows, mixing layers, shocks, and 
thermodynamic and chemical non-equilibrium flow physics. More specifically, scramjet flowpaths are required to 
operate at or near optimum efficiency with minimal aerodynamic loss in flow regimes where flow and chemical time 
scales are of the same order. Consequently, the interaction of turbulence and finite-rate chemical kinetic effects will 
play a leading role in determining flow and flame characteristics and must be modeled with high fidelity. Such 
similarity of scale can also cause the flow regime to be highly unstable and susceptible to flame extinction and re- 
ignition phenomena which can have a detrimental impact on overall flowpath operability and lead to a loss of margin 
and possibly leading to unstart. Unfortunately, due to these multi-scale, multi-physics interactions, this flow regime is 
extremely difficult to model because no assumptions can be accurately postulated about the relationships between time 
and length scales of the multi-physics phenomena. Therefore, the dynamic interactions of multi-physics, multi-scale 
phenomena must be simulated rather than modeled. 

A similar need for the predictive simulation for multi-scale flows in incompressible and subsonic flow has led 
to the development of the large eddy simulation (LES) approach. In LES, the turbulence is separated into large-scale, 
grid-resolved, geometry influenced flow motions, and small-scale, unresolved motions that exhibit some universal 
behavior. 1 This “divide-and-conquer” approach utilizes brute force computing to resolve the most challenging 
dynamical interactions of the flow and relies on “sub-grid” models for the more universal small-scale flows. 
Consequently, LES is inherently 3-D, unsteady and resource intensive. Nevertheless, over the last two decades, 
significant advancements in computational hardware and numerical algorithms have enabled an ability to employ LES 
for the simulation of some reacting flows in practical devices. 2,3 LES also offers a substantial computational cost 
reduction compared with the direct numerical simulation (DNS) approach which has a computational effort that is 
proportional to the cube power of the sub-grid Reynolds number. The value of the sub-grid Reynolds number for a 
typical LES of practical interest varies from 5-100, depending on the characteristic flow Reynolds number, resulting in 
a dramatic two-to-six orders of magnitude reduction in computational cost (when compared with DNS) . 

However, the computational resource requirements of LES are particularly daunting when considering the 
complex, shock dominated, high Reynolds number, and wall-bounded flows found in scramjet flowpaths. The presence 
of strong shocks, shock-shock, shock-turbulence, and shock-boundary layer interactions introduces multi-scale 
interactions that are challenging to resolve numerically. In particular, the number of grid points required to resolve 
near-wall viscous flows via LES is generally prohibitive. The challenges for near-wall viscous flows are similar to 
those for reacting flows because the viscous interactions occur in the sub-grid of LES. For this reason, the majority of 
scramjet LES simulations at NASA Langley Research Center (LaRC) to date have employed a hybridization of 
Reynolds Averaged Navier-Stokes Simulation (RANS) methods and LES methods 4 with advection schemes that utilize 
spatially and temporally 2 nd -order accurate conventional shock capturing upwind methods. 

The computational cost of such a hybrid RANS/LES of a scramjet flowpaths is made up of many 
components. Among these are, the grid density required to spatially resolve the advection and diffusion of the small- 
scale structures, the number of governing equations required to adequately simulate the chemical kinetic processes, 
and the cost of the method used to model the aforementioned physical phenomenon at the sub-grid scale(s). Together 
with reducing the number of transport equations, the reduction of the grid resolution (the total number of grid points, 
cells or elements required to achieve a given level of resolution) is the most direct way to reduce the computational 
resources required. Higher-order (order of accuracy greater than 2) accurate methods have the potential to reduce the 
number of computational cells by an order of magnitude or more. 5 For example, Lele 5 has shown that increasing the 
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formal order of accuracy of a finite-difference differential operator from 2 nd to 6 th -order allows the cell size to be 
increased by a factor of 4, in each computational direction, while maintaining the same level of dissipative error. Of 
particular importance for high-speed flows with shocks, however, are dispersive (phase speed) errors. These errors 
numerically alter the propagation speed of high-wavenumber phenomena thereby introducing decoupling and 
numerical instabilities in addition to significantly altering the physical behavior of the dynamical system. These 
dispersion errors are also significantly reduced by higher-order methods thereby also improving the accuracy of 
simulations. High-order accuracy is also particularly important for LES because turbulence quantities of interest, such 
as scalar variances or reaction fronts, are strong functions of high wavenumber phenomena (the smallest LES-resolved 
scales). Therefore, it is important that the information contained in those wavenumbers be free of numerical artifacts 
that would interfere with their propagation at the correct speed. However when strong discontinuities exist in the flow 
the higher-order accurate finite-difference linear advection schemes typically used, e.g. compact finite-difference, are 
not robust in the vicinity of the discontinuities and the higher-order accurate non-linear discontinuity-capturing 
advection schemes, e.g. WENO, are too dissipative in the smooth parts of the flow. This conundrum has lead to the 
practice of constructing hybrid-advection schemes where a non-dissipative advection scheme is used in the smooth 
parts of the flow and a dissipative discontinuity-capturing advection scheme is used only near the discontinuities. 
There exists a large body of work where various type of higher-order finite-difference advection schemes are 
hybridized and pertinent examples are explored in Refs. 6 and 7. 

However, all is not necessarily lost with respect to 2 nd -order advection schemes. Further utility has been 
demonstrated for 2 nd -order accurate finite-difference methods through the careful hybridization of a dissipative upwind 
advection scheme and a non-dissipative advection scheme 8 which can be shown to be stable. 9 Similar work employing 
2 nd -order finite volume methods utilizing hybridization to significantly reduce their numerical dissipation have also 
been investigated. One such effort involved the hybridization of a non-dissipative 2 nd -order alpha-heuristic skew 
-symmetric centered flux and a dissipative 2 nd -order Weighted Essentially Non-oscillatory (WENO flux). 10 Another 
effort has utilized a 4 th -order divergence form centered reconstruction advection scheme hybridized with 3 rd -order 
upwind biased advection schemes 11 while still another has employed 4 th and 6 th -order least squares stencils and locally 
energy preserving advection schemes hybridized with 2 nd -order upwind biased advection schemes 12 to significantly 
reduce the truncation error of a nominally 2 nd -order cell-centered finite volume advection scheme. These techniques, 
while still formally 2 nd -order accurate on general curvilinear grids, have demonstrated a marked decrease in numerical 
dissipation and a corresponding improvement in resolution of the small scales making such an approach a viable 
alternative worthy of further consideration. Such an approach also offers a way to “upgrade” existing 2 nd -order flow 
solvers such as the VULCAN 13 code developed at the NASA LaRC, thereby reducing numerical dissipation when 
performing hybrid RANS/LES of high-speed propulsion flows. This approach will be explored further in this paper by 
implementing and testing several alternative advection schemes. The schemes considered will include hybridization of 
a 4 th -order non-dissipative advection scheme written as a linear combination of divergence and product rule forms 14 or 
hybridization of a 4 th -order non-dissipative advection scheme based on symmetric reconstruction, 15 with the 
MUSCL, 16 PPM, 17 or 4 th -order WENO 18,19 upwind advection schemes. The hybridization will be accomplished 
employing a “discontinuity sensor” methodology inspired by the approach of Ducros 20 or a sensor proposed by 
Larsson. 21 A series of benchmark problems with increasing spatial and fluid dynamical complexity will be utilized to 
examine the ability of the hybridized advection schemes to resolve and propagate small scale structures typical of 
turbulent flows. The behavior of the discontinuity sensor and it's effect on the detection and capture of strong shocks 
will also be examined. Finally a realistic geometry typically found in a high-speed propulsion system flowpath will be 
computed using one of the examined schemes and the solutions compared with available experimental data. 

II. Numerical Formulation 
A. Smooth flow, non-dissipative advection scheme 

The advection terms in the Euler equations can be written in several alternative forms 14 . Of these forms there 
is significant evidence in the literature that suggests that non-conservative forms that preserve important invariants of 
the continuous equations are preferable to conventional approaches based on the conservative divergence form. 
However, non-conservative forms have historically been viewed as inappropriate for use in flows that contain 
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discontinuities such as those that can arise when the Mach number is greater than one. Fortunately, Fisher et al. 14 
demonstrated that the alternative forms can be recast into an “equivalent, consistent and conservative form thereby 
allowing the Lax-Wendroff theorem to guarantee that the captured discontinuities are weak solutions of the governing 
equations” 14 . These recast alternative flux forms, F a , can be written using a linear combination of the conservative 
divergence form flux, F c , and product rule form fluxes, F e as 

F a = a F c + (1 — a) F e (1) 


where a is the flux combination coefficient. Eq. (1) can be used to write the nonlinear terms in the Euler equations in 
an a — split flux form where the continuity equation is restricted to the divergence form 
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Note that by setting « m =cx e =cx k = cx p = 1 the conservative divergence form of the equations are recovered. 
Other schemes can be obtained using specific values of a that are designed to preserve important invariants of the 
continuous equations. Two examples are the so-called kinetic energy preserving form 12 obtained by setting cx m =l/2 
and <x e = a k = a p = l and the conservative skew-symmetric form of Honein and Moin 22 obtained by setting 
a m =a e =l/2 and « k =cx p =0. Fisher et al. 14 observed that the discrete form of Eqs. (2-4) do not explicitly lead to a 
conservative form unless the following conditions are met: 

1) The discrete spatial operator must be equivalent to a telescoping form. 

2) The discrete spatial operator is Lax-consistent with the physical flux, where the physical flux satisfies the 
Rankine-Hugoniot relations. 

Fisher et al. further demonstrate in Ref. 14 that finite domain ex — split form central-difference operators are 
conservative for any value of ex if a diagonal-norm, finite-domain summation by parts operator is used. They then 
derive discrete 4 th -, 6 th - and 8 th -order summation by parts finite-domain central-difference operators, referred to as the 
2-4-2, 3-6-3, and 4-8-4 operators respectively, that satisfy the aforementioned conditions. The reader is referred to their 
work for the formal proof and stencil coefficients. Due to the control volume heritage of the VULCAN code only the 
interior part of the 4 th - and 6 th -order finite-difference operators were implemented. The interior operators were re- 
interpreted in the cell-center control volume context of the code, i.e. the solution at the cell centers are treated as finite- 
difference solution points and the fluxes at the cell faces are treated as finite-difference flux points. The boundary 


5 

American Institute of Aeronautics and Astronautics 



operators were treated in the same fashion as the interior operators in the manner typical of cell centered finite volume 
codes, i.e. via ghost cells where the number of ghost cells is one half of the stencil width resulting in 2 ghost cells for 
4 th -order and 3 ghost cells for 6 th -order. When a block boundary is a physical wall boundary condition, the ghost cells 
are filled using only the interior data resulting in a “one-sided” stencil near the boundary. When the block boundaries 
connect to other blocks (or domains) the ghost cells are filled by exchanging data with those blocks such that the 
operator is continuous across the block boundaries. This boundary implementation precludes a proof as to the formal 
stability of the approach but it will be demonstrated though numerical testing that the approach has been robust for the 
test cases investigated to date. 

Use of the conservative skew-symmetric a — split form central difference operators has been shown to result 
in a scheme that is stable at very high Reynolds numbers without the introduction of artificial dissipation 23 in the 
absence of discontinuities. However, if discontinuities are present the scheme is not as robust and artificial dissipation 
must be added in the vicinity of the discontinuities or the scheme must be hybridized with a scheme capable of 
capturing discontinuities. There is currently an effort at NASA LaRC to pursue a unified approach where artificial 
dissipation is added only at discontinuities via a finite-difference WENO method. This approach, while being 
mathematically elegant, can be rather expensive and may be difficult to extend to arbitrary mixtures of non-calorically 
perfect gases. Therefore in the current work we elected to hybridize the ex — split flux with the conventional upwind, 
approximate-Riemann-solver-based, discontinuity-capturing methods that already exist in the VULCAN code, i.e. 
fluxes computed using MUSCL, PPM or WENO primitive variable reconstruction with an approximate Riemann 
solver. 


A less expensive alternative non-dissipative approach has also been implemented in the manner 
described by Gieseking et al. 15 where the left and right states are reconstructed to higher-order (in 1-D) at the control 
volume face using a 4 point symmetric interpolation function of the form 


, K L _ K _ -L ^ , 
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q j + I^ q j 
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12 q j+ 2 


( 5 ) 


where q^, denotes the symmetrically reconstructed left, L, and right, R, primitive variables. The non-dissipative 
cell face symmetric reconstruction flux of the conserved quantity Q, F SR (Q), can then be computed in the manner 
typical of compressible flow cell centered finite volume codes from 

FsR(Q)j+l/2 = F (QsR )j+l/2 = ^'(f'( C lsR)j+l/2"^^( C lsR)j+l/2) (6) 


where the stencil notation is illustrated in Fig. 1. A hybrid-advection flux can also be reconstructed using the 
symmetric reconstruction flux, or it can be computed using the symmetric reconstruction states. 

B. Discontinuity capturing, dissipative upwind advection schemes 

The conservative divergence form of Eqs. (2-4), cx m =cx e = (x k = « p = 1, are discretized with a 2 nd -order cell- 
centered, finite volume method. The scheme is nominally 2 nd -order accurate having one quadrature point per cell face. 
A one-dimensional primitive-variable reconstruction at each cell face is performed where the upwind flux can be 
computed using an approximate Riemann solver based approach 

F U p(Q) j+ / 2 =^(F(q E ) j+1/2 +F(q S ) j+1/2 )-^|A(q t ,q S )| j+1/2 (Q(q S ) j+1/2 -Q(q I ) j+1/2 ) ( 7 ) 

where F(q L ) and F(q R ) are the “flux-limited” inviscid fluxes, Q(q L ) and Q(q R ) are the conserved variables 
using the “limited” reconstructed left, L, and right, R, primitive variables, q L and q R , and | A| is computed using 
either Roe's approximate Riemann 24 solver or a Local Lax-Friedricks 25 method using Roe averaged variables. 24 
Alternative upwind advection schemes are also available such as the Low Dissipation Flux Split Scheme 26 or the 
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HLLC flux scheme. 27,29 

Higher-order reconstruction of the primitive variables is accomplished using either a MUSCL, 16 PPM 17 or 
WENO 18,19 scheme. The primitive variables static density, velocity and static pressure are used with the MUSCL 
scheme while the primitive variables static density, velocity and static temperature were found to produce more 
monotone solutions when utilizing the PPM and WENO schemes. When computing mixtures of thermally perfect 
gases the species mass fractions are also used in all three schemes. 

The MUSCL and PPM schemes use gradient limiters during the reconstruction process as the means to 
capture discontinuities. The MUSCL scheme has many limiter options available. However, when utilizing the MUSCL 
scheme as part of a hybrid-advection flux ac = 1/3 scheme with SONIC-A or SONIC-B gradient limiters 28 are 
typically used resulting in a scheme that is nominally 3 rd -order (in 1-D) in smooth flows and lower-order at 
discontinuities. When utilizing the PPM scheme as part of a hybrid-advection flux and if very strong discontinuities are 
anticipated the robustness of the scheme is augmented through the introduction of a pressure based limiter of the form 


j= — max (0, min (l,d 


P j+i-2Pj+Pj-i 

\P j+D +2P j +P j _ 1/ 


+ (<V<5p 


l P j+l~ P jl + l P j~ P j-ll \\ 
min(P j+1 jPpPj-i) " 


( 8 ) 


where P is the static pressure and 5 is a scaling coefficient that typically has a value on the order of 0.1. The 
pressure limiter coefficient, <£ P;j , is computed at the cell centers prior to performing the reconstruction and is used as 
a reconstruction post-processing step that is applied to the PPM reconstructed left and right primitive variables as 

(q L ) j+ D/2=qj+0p,j((q L )j + i/2-qj) (9) 

(q )j+D/2 — q j-t- 1 4 >p, j+ i(q j+ i (q )j+i/2) 

resulting in a scheme that is nominally 4 th -order (in 1-D) in smooth flows and l st -order at discontinuities. 

The four-point WENO reconstruction used in the current work is nominally non-dissipative. However, in the 
presence of a discontinuity, the interpolation of the average states in each control volume (the solution points) to the 
cell interface (the flux point) is modified to ensure that the data is not interpolated across the discontinuity. The stencil 
biasing mechanics of Fisher et al. in Ref. 19 are applied to the primitive variables (q) of the left, center and right 
stencils, S R , S c , and S R respectively shown in Fig. 1. 
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Figure 1. 4 th -order WENO stencil centered about the flux point j+1/2. 

The left state is approximated by, 
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( 11 ) 


3 1 \ 

Qs L — 1 — Qj 2^- 1 ] ’ ^"Qjj > ^ls R ^ ^' c lj+2 _ ^^' c lj+i| 

The weights are computed using 


; z« r 


a =d r 




where d r are the target weights, and the smoothness indicators t and , are defined as 
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(14) 


To prevent the downwind stencil from being assigned a larger weight than the upwind stencils, which would lead to 
instability, the downwind stencil biasing indicator is modified using 




(15) 


The right state at the cell face (flux point) is computed by mirroring the above procedure about the flux point. 

C. Hybridization Methodology and Discontinuity Sensor 

A hybridized cell-face inviscid flux, F H of the conserved quantity, Q, can be computed as a linear 
combination of the upwind flux, and the a — split flux of that quantity based on the stencil discontinuity sensor 
coefficient, (p S s as 

F H (Q) j+ i/D = ^ssF UP (Q)j + i /2 + (l-^ss)F w (Q) j+ i / 2 ( 16 ) 

or as a linear combination of the MUSCL, PPM or WENO reconstructed states, q L,R and the symmetric 
reconstruction states , q^ based on the stencil discontinuity sensor coefficient as 

Qh - <P ss q ~K — ^ssMsr (17) 

When the hybridized states are substituted into Eq. (7), the resulting hybrid symmetric reconstruction flux has the form 

FH(Q)j+i/2 = ^"(F(qH)j+i/2+F (q H )j +1 / 2 ) — ~^\ A (q H , q H )| j+ / 2 ( Q ( qn )j+ 1/2 Q ( Hh) j+ 1/2) (18) 

It is important to note that Eq. (18) will recover the non-dissipative symmetric reconstruction flux of Eq. (6) when the 
discontinuity sensor is zero because q^’ R = q^ thereby causing the dissipation term in Eq. (18) to be identically 
zero recovering F H (Q) j+1/2 =l/2(F(q R ) j+1/2 +F(q R ) j+1/ ) , i.e. Eq. (6). 

The hybrid-advection fluxes are utilized in the VULCAN code within the existing finite volume context 
where the governing equations are written in integral form as 

^T+y / F H (Q)df2 = 0 (19) 

If trapezoidal rule is used to discretize Eq. (19), and the right hand side of Eq. (16) is substituted for F h (Q) in the 
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resulting equation, the discrete form of the equation for the hybrid ex - split flux based approach is obtained as 

[^ssFup(Q)+(l-^ s k s )F^(Q)]dA k = 0 (20) 

where the superscript k denotes evaluation at the k t h face of the control volume. The discrete form of the hybrid 
symmetric-reconstruction based flux approach is obtained in the same manner using Eq. (18) instead of Eq. (16), 
yielding 


SQ 

dt 



-(F k (q^) + F k (qS))--|A k (q^,q^)|(Q(q^)-Q(qJ i )) dA k = 0 


( 21 ) 


Many different approaches exist in the literature for computing discontinuity sensors. Since our intent is to 
use the hybrid-advection scheme in an LES context, we selected a sensor that will cause the use of the non-dissipative 
advection scheme in the vicinity of resolved turbulent content while forcing the use of the dissipative discontinuity- 
capturing advection scheme in the vicinity of shocks. Therefore, in the current work, we chose a discontinuity sensor 
that is a function of the divergence of the velocity vector, V*Uj and the magnitude of the vorticity, | V XuJ. Two 
forms of the sensor were implemented, the Ducros 20 sensor and a sensor due to Larsson 21 . The Ducros sensor has the 
form 


^ss=- 


(V-u ,) 2 

V-Uj) 2 +|Vx Ui |+e D 


( 22 ) 


where e D is small number to prevent division by zero. The Larsson sensor has the form 


where 


^ss - 0 if(0 L <l) 

(// ss =l if (@ L >1) 


0 L =- 


-(V- Ui ) 
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(23) 


(24) 


Here c is the speed of sound, and ds is a cell length scale, e.g. ds = \/ volume cell in 3-D. The constants L ss , i and L ss ,2 
are problem dependent and, for the purposes of the current work, have the values of 5.0 and 0.05, respectively. The 
gradients of the velocity field are computed at the cell centers using the Green-Gauss theorem. Due to the multi-block 
and multi-processor nature of the computational domain, a multi-step sensor post-processing process was adopted 
when computing the Larsson discontinuity sensor where the steps are: 


1) Compute the value of the discontinuity sensor in all the interior cells in the computational domain. 

2) Enforce the discontinuity sensor boundary conditions in the boundary ghost cells. 

3) Communicate the discontinuity sensor between all of the blocks. 

4) Make the value of the discontinuity sensor the maximum of the discontinuity sensor in the 
current cell and nearest neighbor cells. 

5) Enforce the discontinuity sensor boundary conditions. 

6) Communicate the discontinuity sensor between all of the blocks. 

7) Repeat steps 4-6 to increase the domain of influence of the sensor. 

8) Explicitly smooth the discontinuity senor in each block using a Laplacian smoother. 

9) Enforce the discontinuity sensor boundary conditions. 

10) Communicate the discontinuity sensor between all of the blocks. 
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This process ensures that the discontinuity sensor is continuous across block-to-block interfaces, well behaved at 
physical boundaries, and increases the number of cells in the vicinity of shocks that are tagged as shocked cells thereby 
increasing the robustness of the hybrid-advection flux scheme. 

D. Viscous models and discretization 

The form of the Reynolds Averaged Navier-Stokes equations as well as a description of the transport and 
thermodynamic property models employed in the VULCAN code can be found in Ref. 13. The cell face gradients 
required to compute the viscous terms are constructed using the Green-Gauss theorem where the integration is 
performed over an an auxiliary control volume that straddles each cell face resulting in a 2 nd -order accurate 
computation of the viscous fluxes. 

E. Time advancement schemes 

Explicit 3-step 3 rd -order TVD 29 and 5-step 4 th -order 30 Runge-Kutta time advancement schemes as well as 2 nd - 
order backward Euler dual time implicit time advancement schemes are available in the VULCAN code. 

F. Hybrid-Reynolds Averaged Navier-Stokes/Large Eddy Simulation scheme 

The hybrid-Reynolds Averaged Navier-Stokes/Large Eddy Simulation (RANS/LES) approach used in the 
current work is described in detail in Ref. 4. This framework was designed to enforce a Reynolds Averaged behavior 
near solid surfaces that smoothly switches to a LES behavior in the outer portion of the boundary layer and free shear 
regions. Hence, this formulation can be thought of as a wall-modeled LES approach, for which the RANS equations 
are used as the near-wall model. The basic idea is to blend any trusted Reynolds averaged eddy viscosity with a desired 
LES sub-grid scale (SGS) viscosity, along with any transport equations that involve a common RANS and SGS 
property. In this effort, the Menter baseline k — co RANS turbulence model 31 was blended with the one equation SGS 
model of Yoshizawa and Horiuti 32 . 


III. Benchmark Problems 

A series of benchmark problems were investigated to determine the accuracy and examine the behavior of the 
described hybrid-advection schemes. These benchmark problems were divided into five categories: 1) continuous 
inviscid flow, 2) continuous turbulent flow, 3) discontinuous inviscid flow, 4) discontinuous viscous flow containing 
multiple shocks and vortices and 5) discontinuous turbulent flow containing features found in a high-speed propulsion 
flowpath. The continuous inviscid flow problems considered were: the advection of a 2-D isentropic vortex in a Mach 
0.1 freestream 33 and the 3-D inviscid Taylor-Green vortex problem. 34 The continuous turbulent flow problem 
considered was a low Reynolds number, fully-developed turbulent flow in a channel. 35 " 37 These continuous flow 
problems were chosen to test the smooth flow, non-dissipative advection schemes. The discontinuous inviscid flow 
problem was Sod's shock tube problem 38 which was chosen to test the dissipative discontinuity-capturing advection 
schemes. The viscous flow example with multiple shocks and vortices was a viscous version of the inviscid Mach 3.5 
2-D cylinder in crossflow benchmark problem described by Chuadhuri et al. 40 which was chosen to test the 
discontinuity sensors and the hybrid-advection schemes ability to capture complex shock/vortex interactions in a 
temporally evolving viscous flow. The discontinuous turbulent flow containing features found in a high-speed 
propulsion flowpath was supersonic cavity-ramp experiment of Settles et al. 42 

A. Continuous flow benchmarks 

A.1 Advection of a 2-D isentropic vortex 

A grid refinement study of the advection of a 2-D isentropic vortex in a Mach 0.1 calorically perfect, 
y = 1.4, freestream was performed by solving the compressible Euler equations on a square computational domain 
using the finite-volume symmetric-reconstruction-flux scheme (FV-SRF) and the summation-by-parts skew-symmetric 
a — split flux scheme (SBP-ASF). A sequence of three grids, consisting of a 64 x 64 cell coarse grid, a 128 x 128 cell 
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medium grid and a 256 x 256 cell fine grid were used. Each grid contained square cells where the medium grid and 
fine grids had a factor of 2 refinement in each coordinate direction of the coarse and medium grids, respectively. The 
computational domain was solved with periodic boundary conditions. The flow was initialized using the isentropic 
vortex equations of Ref. 33 and the flow was evolved in time for one flow through time of the computational domain 
using a 5 stage 4 th -order Runge-Kutta scheme 30 such that the vortex returned to its initial position on the grid. Figure 2a 
presents fine grid contours of the product of the density and the specific thermodynamic entropy, p S , at time t=0 
where the specific thermodynamic entropy is , 


In (t‘ y" 1 ’ ) (25) 

ln(p) 

and T is the static temperature. The accuracy of each scheme was determined by computing the error, (p S) error at 
every cell center as 

( P S ) error = {P S\ = 0~ (p $)t=A t periodic (26) 

Both schemes provided higher than 2 nd ' order convergence of the error norms with the SBP-ASF being nearly 4 th -order 
on the medium and fine grids. However, the accuracy of the FV-SRF scheme is predominately less than 4 th -order and 
the order property appears to degrade as the grid is refined. 



a) Initial contours 



Figure 2. Fine grid initial contours and profile at y=0 of pS vs. x at time t=0. 


Table 1. Grid refinement study tabulated error analysis for the vortex advection benchmark. 


Size of the Computational Grid 

FV-SRF scheme error 

SBP-ASF scheme error 

L 1 

L 2 

V 

L 1 

L 2 

V 

64x64 

l 

l 

l 

l 

l 

l 

l 

l 

l 

l 

i 

i 

i 

i 

i 

l 

l 

l 

l 

l 

l 

l 

l 

l 

l 

i 

i 

i 

i 

128 x 128 

4.092 

3.778 

3.687 

4.991 

3.879 

4.036 

256x256 

2.753 

3.318 

3.443 

4.121 

3.959 

3.952 


Figure 3 graphically illustrates the error trends of the two low dissipation schemes where a) and b) present the results 
for the FV-SRF and the SBP-ASF schemes, respectively, while Table 1 presents a quantitative measure of the for- 
mal order of accuracy. Further insight can be gained by plotting the error distribution along a y=constant line cut 
through the center of the vortex as denoted by the heavy dashed line in Fig. 2a. Figures 4a, 4c and 4e present the error 
distributions for the FV-SRF scheme on the coarse, medium and fine grids, respectively. Figures 4b), 4d), and 4f) 
present the error distribution for the SBP-ASF scheme on the coarse, medium and fine grids, respectively. The error, 
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represented by the solid red line in these figures, shows that while both schemes have similar levels and tends in the re- 
duction of the error with grid refinement, the level and trends in the reduction of the minimum error outside of the 



a) FV-SRF low dissipation advection scheme accuracy 



a) SBP-ASF low dissipation advection scheme accuracy 


Figure 3. Formal order of accuracy of the non-dissipative schemes for the vortex advection benchmark. 
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10 " 



10 " 


a) Coarse grid, FV-SRF, MUSCL, PPM, and WENO b) Coarse grid, SBP-ASF, MUSCL, PPM, and VVENO 




c) Medium grid, FV-SRF, MUSCL, PPM, and WENO d) Medium grid, SBP-ASF, MUSCL, PPM and WENO 




e) Fine grid, FV-SRF, MUSCL, PPM and WENO 


f) Fine grid, SBP-ASF, MUSCL, PPM and WENO 


Figure 4. Comparison of the p S profiles of all of the advection schemes and the low dissipation scheme error 
on a y=0 line cut through the vortex center on the coarse, medium, and fine grids. 
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vortex is quite different. On the medium and the coarse grid the SBP-ASF scheme error outside of the vortex can 
beseen to be noticeably lower and more uniform than the error of the FV-SRF scheme particularity from x=5 to 
x=10. This indicates that the SBP-ASF scheme should be more accurate for well resolved flow features than the 
FV-SRF scheme. The initial (solid circles) and the final (black lines) p S profiles of all the dissipative and non- 
dissipative advection schemes are also plotted in Fig. 4 and indicate that both schemes exhibit similar phase error 
trends with grid refinement. Figure 4 also highlights how poorly the MUSCL scheme resolves the vortex, relative to 
the other schemes, on all the grids. 

A.2. Inviscid Taylor-Green vortex 

The inviscid Taylor-Green vortex 34 evolves from a resolved initial condition continually stretching and 
producing ever smaller structures analogous to the cascade of structures found in turbulent flow. However, due to the 
inviscid nature of this flow, there is no lower bound on the length scale of these structures because there is no 
dissipation of the structures other than that provided by the advection scheme used. Therefore, the inviscid Taylor- 
Green vortex problem as described and utilized by Johnsen et al. 6 was used to examine the stability of the smooth- 
flow, non-dissipative advection schemes for severely under-resolved scales and to test their ability to temporally 
conserve the mean kinetic energy, (pUjUj), and the mean of the density times the specific thermodynamic entropy, 
(p S). The initial conditions described in Ref. 6 were superimposed on a Mach 0.1, calorically perfect, y = 1.4, 
freestream, and the compressible Euler equations were solved on a periodic 2ttx2ttx2tt computational domain 
consisting of 64 x 64 x 64 cells. The governing equations were evolved in time using the 5 stage 4 th -order Runge-Kutta 
scheme for 20,000 time steps until t*(uref/lref)=18.25, where uref is the freestream velocity, lref is the reference length 
scale, and t is the dimensional time. The temporal evolution of the mean quantities were computed by integrating the 
variables of interest over the computational domain every 50 time steps. The mean of any given variable of interest, 
(g), was computed by summing the product of the cell volume and the cell center value of the variable over the 
computational domain and dividing by the total volume of the computational domain, i.e. 

Ncells Ncells 

<g>= Z vol ig/ Z vol i ( 27 ) 

1=1 1=1 


The initial number and scale of the 3-D vortical structures is illustrated in Fig. 5. After integrating the Euler 



Figure 5. Contours of p UjUj and an iso-surface of p iquj colored by p S , illustrating the scale and 
distribution of the 3-D vortical structures of the inviscid Taylor-Green vortex benchmark case initial flow. 
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a) Temporal evolution of the mean kinetic energy 



b) Temporal evolution of the mean (density*spedfic entropy) 

Figure 6. Comparison of the temporal evolution of (pu i u i )/(pu i u i ) t=0 and (pS)/(pS) t = for the MUSCL, 
PPM, 4 th -order WENO, FV-SRF, and SBP-ASF schemes for the inviscid Taylor-Green vortex benchmark case. 

equations to a non-dimensional time of 18.25 the temporal evolution of (pu^) and (p S) for each scheme was 
examined. Figure 6a) and 6b) present the temporal evolution of (pu i u i )/(pu i u i ) t=0 and (pS)/(pS) t=0 of the non- 
dissipative and dissipative schemes, respectively. Figure 6a) reveals that the MUSCL scheme monotonically depletes 
(pu^) illustrating the dissipative nature of the scheme. The FV-SRF scheme was much less dissipative initially. 
However the FV-SRF scheme began to noticeably depart from 1 at approximately t*(uref/lref)=5, slowly decreasing 
until t*(uref/lref)=10 when the (pu^) begins increasing rapidly. The SBP-ASF scheme, however, is shown to 
preserve the (pu^) with only a very small total increase of approximately 0.2 percent at the end of the simulation. 
Figure 6b) reveals that the (pS) behaves as a mirror image of the (pu^), albeit with a much smaller departure 
from unity. As with the (pu^), the evolution of the (p S) reveals that the MUSCL scheme was significantly more 
dissipative than the other schemes. Figure 6b) also shows that although the temporal evolution of (pS) is different 
for the FV-SRF and SBP-ASF schemes, the magnitude is very nearly the same at the end of the simulation. This 
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decrease in entropy in these simulations clearly violates the second law of thermodynamics and may indicate that as 
the smallest physical scales become unresolvable, the lack of a mechanism to dissipate their energy causes the FV-SRF 
and SBP-ASF advection schemes to violate realizability. 



a) MUSCL + SONIC limiter + HLLC flux 


b) PPM +HLLC flux 


c) 4 th -order WENO + LLF flux 


d) 4%rder WENO + HLLC flux 


Figure 7. Iso-surfaces of p UjUj colored by p S illustrating the scale and distribution of the 3-D vortical 
structures in the inviscid Taylor-Green vortex benchmark case at t= 18.25 for the 
dissipative discontinuity-capturing advection schemes. 
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a) FV-SRF 


Figure 8. Iso-surfaces of p Ujiij colored by p S illustrating the scale and distribution of the 3-D vortical 
structures in the inviscid Taylor-Green vortex benchmark case at t= 18.25 for the non-dissipative advection 

schemes. 


Figures 7a-d) present a visualization of the length scales and their distribution throughout the computational 
domain at t=18.25 for the dissipative, discontinuity-capturing schemes. These figures are ordered by increasing 
accuracy and decreasing dissipation. A comparison of these figures reveals that the scales that are visualized by the 
iso-surface become smaller and more uniformly distributed as the scheme resolution increases and as dissipation 
decreases. Figures 8a-b) present a visualization of the length scales and their distribution throughout the 
computational domain at t=18.25 for the non-dissipative FV-SRF and SBP-ASF schemes. The spatial distribution of 
the small-scale structures computed using the SBP-ASF scheme subjectively appear to be visibly more uniform than 
the structures computed using the FV-SRF scheme. Moreover, the spatial variation of p S in the SBP-ASF simulation 
appears to be significantly different and more oscillatory than the FV-SRF simulation. The comparison of Figs. 7 and 8 
serves to illustrate how much more dissipative the upwind schemes, MUSCL, PPM and WENO, are in comparison to 
the low dissipation FV-SRF and SBP-ASF schemes. 


A.3 Fully developed turbulent channel flow, unit Re=5,600.0 


A fully developed turbulent channel flow 35,36 was chosen to test the behavior of each numerical scheme for a 
wall-bounded LES application. This configuration is an ideal wall-bounded flow configuration for testing numerical 
algorithms for LES, since periodic boundary conditions can be applied in both the streamwise and spanwise directions 
(avoiding the uncertainties associated with the specification of inflow/outflow conditions with resolved turbulent 
content). The use of periodic conditions in the streamwise direction requires that the mean pressure gradient (which 
physically drives the flow) be replaced with an equivalent spatially uniform source term that balances the wall shear 
stress at each temporal integration step. A snapshot of the turbulent flowfield is presented in Fig. 9. This image 
displays the resolved turbulent structures along the lower half of the channel using iso-surfaces of vorticity colored by 
Mach number. The boundary planes are contours of vorticity magnitude. Several sources of data at a Reynolds number 
of 5600 (based on the bulk velocity and channel height) are available for comparisons with the present numerical 
results which include measurements (Eckelmann 37 ), DNS (Kim et al. 35 ), and LES (Okong'o and Knight 36 ) data. 
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Figure 9. A snapshot of the turbulent flowfield in a channel displaying the resolved turbulent structures along 
the lower half of the channel using iso-surfaces of vorticity colored by Mach number with the boundary planes 

being displayed using contours of vorticity magnitude. 

The simulations were carried out on a domain of length 2 tt h in the streamwise direction and tt h in the 
spanwise direction (where h represents the channel height). The physical dimensions of this domain were chosen to 
match the DNS of Kim et al. 35 where two-point correlations were analyzed to ensure that the computational domain 
was large enough such that the turbulent fluctuations were uncorrelated at a separation of half the domain width in 
each homogeneous direction. The Mach number based on the bulk velocity was set to 0.25. This value was chosen to 
balance computational efficiency of the compressible solver with the desire to simulate an incompressible flowfield. 
The resulting mean static density variation across the channel was approximately 1%. Hence, comparisons with the 
incompressible DNS data source are meaningful. The closure for the Sub-Grid- Scale (SGS) terms was achieved using 
a constant coefficient Smagorinsky model with a van Driest damping factor to facilitate integration through the inner 
portion of the boundary layer: 


n= f vd C s p^V(2S ij S ij ) 

where 


S = — 

1J 2 


dm dm 
— -+ — 1 
dxj dx i 


. fvd=l-e l 


,(-y+26) 


, y-y(u T /v) 


(28) 


(29) 


In the above expression, f V d is the van Driest wall damping term, C s is the Smagorinsky constant, p is the static 
density, A g is the filter width, and m are the components of the Cartesian velocity vector. The friction velocity, u T , 
required by the van Driest damping term was taken to be the spatially-averaged value over both channel surfaces and 
the filter width, A g was taken to be the local wall-normal cell height, and C s was set to 0.01. All channel surfaces 
were assumed to be isothermal with a wall temperature of 298.15 K. 

The simulations were performed on three different grids with streamwise (x), wall-normal (y), and spanwise 
(z) cell counts of 32 x 64 x 48 (coarse), 64 x 64 x 96 (medium), and 128 x 64 x 96 (fine). The wall-normal distribution 
for each grid had a minimum spacing at the wall of 0.00278 h (Ay ^0.05) , and a maximum spacing of 0.0311 h 
(A z^6) , at the channel centerline. The grid was uniformly spaced in each of the homogeneous directions. The 
finest streamwise grid spacing was Z\x^9, and the finest spanwise grid spacing was Z\z^ . The finest grid 
employed closely matched the grid used in the DNS of Kim et al. 35 which utilized a spectral method (Fourier series in 
the homogeneous directions, and a Chebychev polynomial expansion in the wall-normal direction) to evaluate the 
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spatial derivatives. The medium grid was a close match to the grid used in the LES of Okong'o and Knight 36 . Their 
work employed a finite-volume methodology with a third-order least squares variable reconstruction for the fluxes at 
the cell faces. This spacing was not altered for the various grids considered. As a result, the differences in the statistics 
obtained from each grid is purely a function of the level of resolved turbulence captured by the numerical scheme. 

All simulations were advanced in time using the 5-stage 4 th -order Runge-Kutta scheme with a time-step of 20 
microseconds. The initialization procedure was performed on the coarse grid (via the FV-SRF scheme) using a fully- 
developed laminar profile with the velocity field randomly perturbed as the initial condition. After 1 to 2 eddy-turnover 
times (defined as 0.5 h/u T ) , the perturbed laminar flow begins to transition towards a turbulent state. The transition 
process completed after approximately 5-6 eddy-turnover times, and the solution at this point was saved and used to 
initialize all subsequent simulations. The transition process is illustrated by the temporal evolution of the surface- 
averaged friction velocity in Fig. 10. Starting from the fully transitioned turbulent flowfield obtained through this 
initialization process, five different inviscid flux schemes were used to simulate this flowfield on the 64 x 64 x 96 grid. 
The schemes tested were SBP-ASF, FV-SRF, WENO-HLLC, PPM-HLLC, and MUSCL-HLLC with the SONIC-A 
flux limiter. The resulting surface-averaged friction velocity time histories are shown in Fig. 10. On this grid, the 
dissipation associated with the MUSLC-HLLC scheme was large enough to completely damp the resolved turbulent 
fluctuations; eventually resulting in a nominally laminar profile. The PPM-HLLC and WENO-HLLC schemes were 
able to maintain some of the resolved turbulent content, but clearly a large fraction of this content was dissipated 
through numerical diffusion. Although not shown, the WENO-HLLC scheme was also able to marginally maintain 
some of the resolved turbulent content when the 32 x 64 x 48 grid was used, but the PPM-HLLC scheme was unable to 
maintain any of the resolved turbulence. Both of the non-dissipative schemes were able to maintain levels of resolved 
turbulence comparable to the DNS values from Ref. 35 (even on the coarsest grid used in this effort). Hence, only the 
non-dissipative schemes (SBP-ASF and FV-SRF) were considered for further evaluation. 



Figure 10. The temporal evolution of the surface-averaged friction velocity using the SBP-ASF, FV-SRF, 
WENO-HLLC, PPM-HLLC, and MUSCL-HLLC with the SONIC-A flux limiter advection schemes 
on the 64 x 64 x 96 grid (t_ref is the eddy-turnover time). 

The effect of grid resolution on the surface-averaged friction velocity time histories is shown in Fig. 11. For 
brevity, only the results from the SBP-ASF scheme are shown since the FV-SRF results were similar. The results from 
the 32 x 64 x 48 and 64 x 64 x 96 grids produced somewhat higher friction velocities than the DNS of Kim et al. 
However, the friction velocity obtained from the use of the 128 x 64 x 96 grid compares quite well with the DNS 
values. Prior to the gathering of flowfield statistics, each simulation was allowed to progress an additional 6 eddy- 
turnover times as indicated by the vertical dashed line in Fig. 11. Statistics were then gathered for approximately 50 
eddy turn-over times. The statistical error was further reduced by averaging over each of the spatially homogeneous 
directions to produce averaged properties that are a function of the wall-normal coordinate only. 
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Figure 11. The temporal evolution of the surface-averaged friction velocity using the SBP-ASF and FV-SRF 
advection schemes on the 32 x 64 x 48 (coarse), 64 x 64 x 96 (medium) and 
128 x 64 x 96 (fine) grids (t_ref is the eddy-turnover time). 

The mean velocity profile and turbulence statistics extracted from the 128 x 64 x 96 grid simulations are compared to 
measurements and other independent computational efforts in Figs. 12 and 13. Figure 12 presents the FV-SRF and 
SBP-ASF mean velocity profiles showing that they match the DNS data quite well, but there is a discrepancy noted 
with respect to the measurements in the log layer. Kim et al. 35 noted that these measurements were inconsistent with 
earlier measurements taken at the same facility but at a slightly larger Reynolds number. Their analysis showed that a 
multiplicative "correction factor" of 1.06 on the measured friction velocity was required to match the earlier 
measurements, and this same correction factor matched their DNS data as well. The normalized rms velocity 
fluctuations and Reynolds shear stress are compared in Figs. 13 and 14 where Fig. 13 presents the results of the FV- 
SRF and SBP-ASF advection schemes on the fine grid, and Fig. 14 presents the results from the FV-SRF advection 
scheme on the coarse, medium and fine grids. The level of agreement with the DNS data is quite good for the 
Reynolds shear stress and the rms streamwise velocity fluctuation profiles. When the medium and fine grids are used, 
the current fine grid LES results consistently underpredict the DNS values, but the worst case percent difference is less 
than 5%. Moreover, the current medium grid LES results show a noticeable improvement over the previous LES effort 
performed on grid similar to the current medium grid. 



Figure 12. A comparion of the mean streamwise velocity profile (in wall units) computed using 
the FV-SRF and SBP-ASF advection schemes on the 128 x 64 x 96 grid with 
pre-existing computational and experimental results. 
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y / h y / h 






Figure 13. Comparison of normalized streamwise fluctuating rms velocity (upper left), wall-normal fluctuating 
rms velocity (upper right), spanwise fluctuating rms velocity (lower left), and Reynolds shear stress 
(lower right) computed using the FV-SRF and SBP-ASF advection schemes on the 
128 x 64 x 96 (fine) grid with pre-existing computational results. 
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Figure 14. Comparison of normalized streamwise fluctuating rms velocity (upper left), wall-normal fluctuating 
rms velocity (upper right), spanwise fluctuating rms velocity (lower left), and Reynolds shear stress 
(lower right) computed using the FV-SRF advection scheme on the 32 x 64 x 48 (coarse), 

64 x 64 x 96 (medium), and 128 x 64 x 96 (fine) grids with pre-existing computational results. 
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B. Discontinuous flow benchmarks 
B.l Sods 1-D shock tube 

The 1-D shock tube problem of Sod 38 was computed to examine the behavior of the dissipative discontinuity- 
capturing advection schemes described in section II. B. Figure 15 presents the computational domain and the flow 
initialization. The computational domain was discretized with 100 cells in the tube x-direction and 5 cells in the tube 
y-direction. The computational domain was sub-divided into 2 blocks of 50 x 5 cells each initially separated by an 
impermeable diaphragm, indicated by the dashed line in Fig. 15. Table 2 presents the flow primitive variable 
initialization values for each block. 



block 1 

diaphragm block 2 





x ^ 

Figure 15. Computational domain for Sod's shock tube benchmark case. 


Table 2. Initial flow conditions for Sod's shock tube benchmark case. 


flow primitive variables 

block 1 

block 2 

static density, p , kg/m 

9.47151 

1.18394 

velocity, u, m/sec 

0.0 

0.0 

static pressure, p, Pascals 

1013250.0 

101325.0 


The flow of a calorically perfect, y = 1.4, gas was evolved in time from the initial conditions at t=0, 
assuming that the diaphragm burst instantaneously, to a non-dimensional time of t=0.00721 using a 5 stage 4 th -order 
Runge-Kutta scheme. The MUSCL and PPM reconstruction schemes were run with the HLLC approximate Riemann 
solver, and the WENO reconstruction scheme was run with the HLLC and the local Lax-Friedricks (LLF) 
approximate Riemann solvers to examine the individual effects of the different reconstruction schemes and 
approximate Riemann solvers when capturing the discontinuous features of the flow. Figures 16a-d) present the 
computed and analytic non-dimensional density profiles for each reconstruction/approximate Riemann solver 
combination with Fig. 16a) also identifying the key physical phenomena of the shock tube flow. 

A cursory examination of the shock tube solutions reveals that all three reconstruction schemes give similar 
results for a given approximate Riemann solver. Although not shown, the most monotone behavior for the PPM and 
the 4 th -order WENO schemes in the vicinity of the contact discontinuity was achieved using density, velocity, and 
static temperature as the reconstruction variables while density, velocity and static pressure worked equally well for the 
MUSCL scheme. A more detailed examination reveals that the biggest differences between solutions occur in the 
vicinity of the trailing edge of the rarefaction wave, the contact discontinuity and the shock. Differences at the trailing 
edge of the contact discontinuity appear to be primarily due to the choice of Riemann solver whereas differences at the 
contact discontinuity and the shock appear to be primarily due to the order of accuracy of the reconstruction scheme. 
Overall the PPM and the WENO schemes seem to give slightly better resolution at the trailing edge of the rarefaction 
wave, the contact discontinuity, and the shock than the MUSCL scheme, and the PPM scheme is slightly better than 
the 4 th -order WENO scheme at the trailing edge of the rarefaction wave and the shock. Consequently a choice between 
PPM and the 4 th -order WENO scheme was deferred until more is known about the cost and robustness of the 
individual schemes for a wide variety of problems. It is also important to note that the WENO scheme can and has 
been extended to 6 th -order 39 which, when combined with the 6 th -order SBP-ASF scheme, would enable the 
implementation of a 6 th -order hybrid-advection scheme . Therefore, the PPM and the WENO schemes are preferable to 
the MUSCL scheme because of the slightly improved resolution for discontinuous flow and because their baseline 
order of accuracy, i.e. accuracy when the flow is smooth, are consistent with the nominal 4 th order accuracy of the non- 
dissipative advection schemes. Furthermore, although the HLLC approximate Riemann solver was found to b 
preferable to the LLF approximate Riemann solver, LLF continued to be considered because it is commonly used with 
WENO schemes. 
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Density Density 



a) MUSCL + SONIC limiter + HLLC flux 





c) ^-order WENO + LLF flux d) ^-order WENO + HLLC flux 

Figure 16. Comparison of the computed and analytic non-dimensional density profiles for Sod's shock tube 
benchmark using the MUSCL, PPM and WENO schemes with the 
HLLC and the LLF approximate Riemann solvers. 
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C. Hybrid advection scheme test case 

C.l 2-D supersonic viscous flow past a circular cylinder 

The 2-D Mach 3.5 inviscid flow of a calorically perfect, y = 1.4, gas over a circular cylinder with the flow 
conditions of Chaudhuri et al. from Ref. 40 was simulated to test the ability of selected hybrid-advection schemes to 
resolve a flow rich in shock-shock and shock- vortex interactions. These types of interactions are also of particular 
interest when computing hypersonic propulsion flow paths. To that end, a viscous flow over the Ref. 37 geometry at 
similar conditions and a moderate Reynolds number was computed to provide a complex temporally evolving flow 
that would test the ability of the hybrid advection scheme to robustly and reliably capture and resolve shock-shock and 
shock-vortex interactions. The gas properties and inflow conditions are presented in Table 3. Figure 17 presents 
the computational domain 


Table 3. Gas properties and inflow conditions for the supersonic viscous flow over a circular cylinder 


Gas 

Constant 

y 

Mach no. 

Static 

temperature 

Static 

pressure 

Unit 

Reynolds no. 

Molecular 
Prandtl no. 

287.0 J/(kg-K) 

1.4 

3.5 

298.15 K 

101324.13 Pa 

10,000 

0.70 


and the density contours of a laminar, quasi-steady-state solution that was used as the initial condition for the time 
accurate simulations that were performed to test the hybrid advection scheme options. The boundaries at x/D=-3 and 
x/D= 21 of the computational domain were treated as supersonic inflow and outflow boundaries respectively, the 
boundaries at y/D=+3 and y/D=-3 were treated as inviscid walls, and the cylinder wall was treated as a no-slip 
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Figure 17. The computational domain and density contours of the quasi-steady state solution 
of a Mach 3.5 viscous 2-D flow over a cylinder that was used as the initial flow solution 
for testing the hybrid advection scheme and discontinuity sensors. 
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adiabatic wall. The baseline grid consisted of a total of 327,680 cells that was decomposed using VULCAN's block 
splitting utility into 80 equal sized blocks for efficient parallel processing. The initial solution was obtained using a 
diagonalized approximate factorization scheme with local time stepping at a CFL of 4.5. The PPM reconstruction 
scheme with the HLLC approximate Riemann solver was used to compute the inviscid fluxes and an auxiliary cell 
volume Green-Gauss methodology was used to compute the cell face gradients for the viscous fluxes of the 
compressible full Navier-Stokes form of the governing equations. The time accurate computations were performed 
using the 5 stage 4 th -order Runge-Kutta scheme. The behavior of the discontinuity sensor functions was examined by 
evaluating the functions at the pseudo-steady state initial flow conditions and plotting their contours. The results 
obtained using the Ducros discontinuity sensor function of Eq. (22) and the Larsson discontinuity sensor function of 
Eqs. (23 and 24) are presented in Figs. 18a) and 18b) respectively. 
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Larsson Shock Sensor Contours 
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Figure 18. Discontinuity sensor contours computed using the Ducros and Larsson functions for the the pseudo- 
steady state initial flow solution where blue indicates where the non-dissipative scheme would be used and red 
indicates where the dissipative discontinuity-capturing advection scheme would be used. 
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In addition, the discontinuity sensor multi-step post-processing process described in section II.C was used in 
conjunction with the Larsson sensor. The contour plots in Figs. 18a) and 18b) are interpreted as follows: the areas of 
pure blue are where the sensor has a value of 0, and the areas of pure red are where the sensor has a value of 1. Colors 
in between indicate values between 0 and 1. Inserting these values into Eqs. (16) or (17) reveals that the non- 
dissipative advection scheme would be used in the blue areas and the dissipative discontinuity-capturing upwind 
advection scheme would be used in the red areas and a linear combination of the schemes depending on the 
discontinuity sensor value would be used everywhere else. Examination of Fig. 18a) in this context immediately 
reveals that the Ducros discontinuity sensor is indicating that the non-dissipative advection scheme would be used only 
upstream of the bow shock, in the cylinder boundary layer and in the wake. Examination of Fig. 18b) reveals that the 
Larsson discontinuity sensor is indicating that the non-dissipation advection scheme would be used nearly everywhere 
except at the shocks which appear to be clearly delineated by the sensor. This marked difference in behavior of the 
sensors is due to two reasons: 1) the threshold function in the denominator of the Eq. (24), i.e. the L ss , 2 (c/ds) term, 
prevents the function from activating due to small oscillations in the flow solution and 2) the negative sign in the 
denominator of Eq. (24) and the lower limit function of Eq. (25) forces the discontinuity sensor value to be 0 for 
expanding flows whereas the Ducros discontinuity sensor is equally active in expanding and compressing flow. In 
addition, Eq. (25) acts as a Heaviside step function causing the Larsson discontinuity sensor to have a binary behavior, 
i.e. it only produces values of 0 or 1. However, it is important to reiterate that, in an effort to improve robustness, this 
binary behavior is deliberately smeared by the use of the discontinuity sensor multi-step post-processing approach that 
adds a “fringe” of values between 0 and 1 in the vicinity of the shocks. Overall, the behavior of the Larsson 
discontinuity sensor was found to be preferable for three reasons: 1) the non-linearity of the hybridization process is 
restricted to a small number of cells, 2) the 4 th -order non-dissipative scheme is used more, thereby improving the local 
accuracy, and 3) the expense associated with the dissipative discontinuity-capturing advection scheme is restricted to a 
small number of cells, thereby significantly reducing the total computational effort required to compute the inviscid 
fluxes. 

Based on the results of the preceding benchmark tests, the time evolution of supersonic flow over the cylinder 
was computed using three different hybrid-advection scheme configurations as described in Table 4. These 
solutions were computed in a two step manner. The first step was to evolve the solution in time from the pseudo-steady 
state initial condition for two flow through times of the computational domain to establish temporally periodic flow. 
The second step was to evolve the solution in time approximately eight flow through times while performing a density 
weighted time average, i.e. a Favre average, of the flow solution with a constant sampling frequency. In addition, 
instantaneous flow solutions were gathered with a constant sampling frequency to allow animation of preselected flow 
variables to be constructed. 


Table 4. Hybrid-advection scheme configurations used to compute the cylinder in crossflow case 


Case No. 

Dissipative discontinuity-capturing advection scheme 

Non-dissipative 

advection 

scheme 

Discontinuity 

Sensor 

Higher-order reconstruction scheme 

Approximate Riemann solver 

1 

PPM 

HLLC 

FV-SRF 

Larsson 

2 

PPM 

HLLC 

SBP-ASF 

Larsson 

3 

WENO 

LLF 

SBP-ASF 

Larsson 


Figure 19 presents a sample instantaneous flow solution from case 2 in Table 4 where Fig. 19a presents the 
instantaneous density contours, and Fig. 19b presents the instantaneous discontinuity sensor contours. The density 
contours illustrate the complexity of the shock-shock, shock-shear layer, and shock vortex interactions that occur. The 
discontinuity sensor contours demonstrate that the Larsson discontinuity sensor function has done a reasonably good 
job of detecting the predominant flow discontinuities. However some of the discontinuity sensor contours are 
fragmented, i.e. not continuously delineated. This behavior could possibly result in robustness issues in more 
challenging computations. Therefore finding a better discontinuity sensor is a subject that deserves more attention. 
Johnsen et al. 6 make a similar observation and recommend alternatives approaches that deserve further investigation. 
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b) Instantaneous contours of the Larsson shock sensor 

Figure 19. Behavior of Larsson discontinuity sensor for Mach 3.5 viscous 2-D flow over a cylinder 
using the case 2 (Table 4) hybrid-advection scheme configuration from table. 
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The Favre averaged solutions were collected for the three hybrid-advection scheme solutions and were 
compared with the results from a grid refinement study performed by a High-Order L 2 stable multi-domain finite- 
difference method named HOFD 41 , currently under development at the NASA Langley Research Center for boundary 
layer stability/transition simulation. In that study the HOFD code was used to simulate the supersonic viscous flow 
over the cylinder on a sequence of grids using a 4 th -order WENO method. The grid used in the current work is the 
medium grid from that grid sequence. The coarse and fine grids of the sequence were uniform factor of 2 coarsening 
and refinement of the medium grid, in each coordinate direction, resulting in coarse, medium, and fine grids that 
consisted of 81,920 (coarse grid), 327,680 (medium grid) , and 1,310,720 (fine grid) cells, respectively. Figure 20 
compares the medium grid VULCAN cases 1, 2 and 3 non-dimensional Favre averaged density, p! Pm , solution 
profiles with the HOFD coarse, medium and fine grid solutions along Y/D=0 and X/D=5 cuts through the cylinder 
wake. Figure 20 shows that there are small differences between the VULCAN medium grid solutions and the 
HOFD coarse, medium and fine grid solutions near x/D=5, with the VULCAN medium grid solutions 



X/D 


Figure 20. Comparison of VULCAN cases 1-3 (Table 4) and HOFD Favre averaged p/p re , profiles 
at the Y/D=0 station for Mach 3.5 viscous 2-D flow over a cylinder. 
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Figure 21. Comparison of VULCAN cases 1-3 (Table 4) and HOFD Favre averaged pi P re f profiles 
at the X/D=5 station for Mach 3.5 viscous 2-D flow over a cylinder. 

(4 th -order WENO) being arguably closer to the HOFD fine grid solution than the HOFD medium grid solutions. Figure 
21 shows similar good overall agreement while also revealing a small amount of asymmetry in all the solutions 
indicating that a statistical steady state was not quite reached in the simulations. It should also be noted that the 
interior advection and diffusion schemes implemented in the HOFD code are both 4 th -order accurate, the boundary 
treatments are 3 rd -order accurate, and the grids used were constructed such that the grid in each block was analytic, 
thereby making it possible for the HOFD code to be formally 4 th -order accurate for this grid. In addition, it is also 
notable that HOFD has significantly less dissipation than “classical” WENO schemes due to the high-order L 2 stable 
formulation used. Therefore, the favorable comparison shown in Figs. 20 and 21 is reassuring in that it appears that the 
hybrid-advection scheme is doing a reasonable job of capturing the discontinuities while not introducing excessive 
amounts of numerical dissipation. 
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D. Hypersonic Propulsion Flowpath Application: 

D.l Settles Supersonic Ramp Cavity Hybrid Reynolds Averaged Navier-Stokes Simulation/Large Eddy 
Simulation 

The next test problem used to evaluate the candidate numerical schemes involves a supersonic flow detaching 
from a rearward facing step and reattaching along an inclined wall further downstream. This configuration contains 
many of the physical flow features present in hypersonic propulsion flowpaths (e.g. massive flow separation, and shock 
wave / turbulent boundary layer interactions). Settles et al. 42 conducted experiments using the test article shown in 
Fig. 22. The experiments were performed in the Princeton University 20cm by 20cm High Reynolds Number 



Figure 22. Experimental test article with one aerodynamic fence removed 
and the computational domain outline of Fig. 23 signified by dashed lines. 

Supersonic Wind Tunnel and consisted of Mach 2.92 air entering the tunnel test section at a stagnation temperature of 
258K, a stagnation pressure of 0.69 MPa, and a free stream unit Reynolds number of 6.7xl07/m. The test article was 
designed to allow a turbulent boundary layer to develop along a flat plate before separating over a sharp corner. The 
free shear layer that is formed at this separation point eventually reattaches along a wall of length 18 cm inclined at 20 
degrees. The position of the inclined wall was adjusted for each experimental test run, before taking measurements, to 
minimize changes in cavity pressure (and therefore shear layer deflection) in order to match the pressure of the flow in 
the cavity with that of the freestream. The cavity depth was 2.54 cm, and the model was sized and centered in the test 
section to ensure that it was not directly influenced by the tunnel sidewall boundary layers. In order to further promote 
two-dimensionality, aerodynamic fences lined the edges of the inclined wall. Mean flow measurements were made 
using a traversing set of Pitot and static pressure probes, and these measurements, in conjunction with an assumption 
of constant total temperature, enabled the researchers to calculate Mach number, velocity, density, and other relevant 
flow properties. A spark shadowgraph enabled observation of major flow features. 

This experiment was recently simulated by Quinlan et al. 43 using the aforementioned VULCAN hybrid 
LES/RANS framework with the FV-SRF/PPM-HLLC hybrid advection scheme and the Ducros and Larsson 
discontinuity sensors. Table 5 presents a summary of the computational results performed in Ref. 43 that compared 
select hybrid-advection schemes described in current work with the experimental data of Settles 42 . 


Table 5. A summary of the four simulations extracted from the work of Quinlan et al. 42 and reported herein. 


Case Name 

Method 

Inflow Treatment 

No. Grid Cells 

Advection Scheme 

Discontinuity Sensor 

R35D 

LES/RANS 

Recycling 

35.0 x 10 6 

FV-SRF/PPM-HLLC 

Ducros 

R35L 

LES/RANS 

Recycling 

35.0 x 10 6 

FV-SRF/PPM-HLLC 

Larsson 

R8D 

LES/RANS 

Recycling 

8.8 x 10 6 

FV-SRF/PPM-HLLC 

Ducros 

NR8RAS 

RANS 

No Recycling 

8.8 x 10 6 

PPM-HLLC 

N/A 
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In order to introduce coherent turbulent structures at the inflow, turbulent fluctuations and relevant turbulence 
properties are recycled and rescaled from a plane one cavity-depth upstream of the step. The turbulent fluctuations are 
extracted from the donor plane by subtracting time- and span-averaged mean flow data from the instantaneous flow 
data. These fluctuations are scaled according to boundary-layer similarity laws before being superimposed to a mean 
flow profile. The general procedure used for this process is described by Choi et al. 45 . 

Simulations were performed on a structured, three-dimensional grid consisting of approximately 35 million 
cells. The grid was composed of four primary blocks, as shown in Fig. 23. The two upper blocks are connected to the 



coarsened once in each coordinate direction for visual clarity. 

lower blocks via non-CO continuous patches, allowing the farfield boundary to be placed far from the domain of 
interest without significantly effecting the computational costs. The lower upstream block is attached to the lower 
downstream block by a simple non-CO continuous patch where only the spanwise coordinates are mismatched (the 
downstream spacing is twice that of the upstream spacing). This strategy was based on the premise that the grid 
resolution requirements in the spanwise direction would not be as stringent as that for the attached wall boundary layer. 
The streamwise, transverse and spanwise dimensions of the blocks adjacent to the wall directly upstream and 
downstream of the step were 525 x 121 x 301 and 453 x 241 x 151, respectively. The dimensions of the upstream and 
downstream farfield blocks were 29 x 45 x 15 and 153 x 45 x 15, respectively. The patches connecting the lower and 
upper blocks were positioned far enough from the domain of interest to ensure that resolved boundary layer and free 
shear layer eddies are not convected into them. The physical dimensions in the x, y, and z directions of the 
computational domain are 28.6cm, 22.5cm, and 3.81cm, respectively. The flat plate extended four step heights 
upstream of the step to allow sufficient space for recycling of turbulent fluctuations, and the grid was clustered to 
provide a nominal y+ of 1.0 adjacent to the surface. The streamwise and spanwise grid spacing was set to l/20th of the 
nominal approach flow boundary layer thickness. Downstream of the step, the grid clustering was relaxed toward a 
more isotropic grid for the resulting free shear layer. The maximum x, y, and z spacing in this region was 2.5%, 5.0%, 
and 1.0%, respectively, of the cavity depth (2.54 cm). The cell aspect ratio for the LES regions of the domain was no 
greater than approximately 20.0. The three-dimensional grid was split into 1980 sub-blocks to allow efficient run-time 
parallelization on the NASA Advanced Supercomputing resources. Each simulation required approximately 240 hours 
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of wall clock time on 512 processors. 

Mean inflow variable profiles and flowfield initializations for the simulations were obtained using a 
converged RANS simulation. The upstream domain for the RANS solution was sized to allow the freestream flow to 
develop until the momentum thickness one cavity-depth upstream of the step matched the experimental measurements. 
The nominal freamstream had a static pressure of 21.2 kPa, Mach number of 2.92, and temperature of 92.4 K. The 
hybrid LES/RANS simulations utilized periodic boundary conditions at the z-min and z-max boundaries. All walls 
were assumed adiabatic and no-slip, an extrapolation condition was applied at the exit plane, and a characteristic 
condition was enforced at the far-field plane. In order to expedite the formation of resolved turbulent content, initial 
turbulent fluctuations were superimposed onto the initial Reynolds-averaged solution upstream of the cavity step. 
These turbulent fluctuations were extracted from a previously-simulated (via LES) generic turbulent boundary layer on 
a flat plate and rescaled to match the integral characteristics of the approach boundary layer. 

Based on the results from the previous benchmark problems, the hybrid-advection schemes chosen were the 
PPM+HLLC/FV-SRF and the PPM+HLLC/SBP-ASF schemes. However difficulties were encountered running the 
PPM+HLLC/SBP-ASF hybrid-advection scheme. The cause of these difficulties was subsequently determined but the 
computations could not be completed in time to be included in the current work. The switching function forms of 
Ducros and Larsson were both considered in this effort (see Table 5 for the case identifiers). The governing equations 
were integrated in time using a dual-time-stepping approach, in which a diagonalized approximate factorization (DAF) 
method was used for integration in pseudo-time with a maximum Courant-Friedrichs-Lewy (CFL) number of 100.0, 
and a second-order three-point backwards finite-difference approximation was used for integration in physical time. 
Sub-iteration convergence was considered satisfied when a residual reduction of 2.5 orders-of-magnitude was met or a 
limiting number of 15 sub-iterations. The physical time step was set to 0.1 ps. Following initialization, the flowfield 
was evolved for a minimum of 10 characteristic flow-through times before taking statistics. A characteristic flow- 
through time was defined as the time required for a particle in the free stream to travel from the corner of the step to 
the point of reattachment along the inclined wall (approximately 0.2ms). Once a statistically stationary state was 
achieved, flowfield properties were ensemble-averaged in time at a rate of 0.5 microsecond for periods of 5, 10, and 
15 flow-through times and subsequently averaged spatially in the homogeneous spanwise direction. 

Prior to comparing results to experimental data, the approach flow wall boundary layer was compared to law- 
of-the-wall theory to confirm that the wall boundary layer developed properly upstream of the step. This confirmation 
was necessary to ensure that a constant required by a blending function used by the hybrid LES/RANS model was 
chosen appropriately, and that the recycling procedure itself did not adversely affect the boundary layer development. 
Boundary layer velocity profiles at several planes upstream of the step were extracted and compared to theoretical 
predictions and the steady-state RANS solutions using the van Driest II transformation. The extracted planes included 
the inflow plane, the donor recycling plane, and the plane at the step corner. For each plane, the averaged LES/RANS 



Figure 23. Streamwise velocity profile (in wall units) and the blending function profile 
at the recycling plane, shown in Fig. 24, of the cases in Table 5. 
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results agreed well with both the theoretical log-layer and with steady-state RANS data. The comparison of the profiles 
at the recycling plane are included in Fig. 23 to demonstrate the level of agreement. At each plane, the LES/RANS 
simulations reproduced the slope and location of the log-layer well, and the blending function anchored the mean 
LES/RANS transition point near the upper edge of the log-layer as intended. 

Figure 24 provides several snapshots of instantaneous velocity contours to show the major flow features 
anticipated for this flowfield. The supersonic flow upstream of the cavity detaches at the step corner, thereby creating a 
free shear layer over the cavity. This free shear layer attaches along the inclined wall and interacts with an oblique 
shock front standing off the inclined wall. Fluid is entrained into the cavity near the location of shear layer attachment, 
and this fluid forms a recirculation zone which drives an unsteady low-frequency motion of the free shear layer. Note 
that annotations are included to indicate major flow features and identify each image's relative location in time, which 
is given as the value of elapsed time, A t , relative to an arbitrary origin, divided by the characteristic flow-through 
time, A t . Also indicated are the locations of the measurement planes where comparisons will be made to 
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Figure 24. Instantaneous snapshots of velocity magnitude at 5 time slices for case R35D (Table 5) 
(recycle plane, major flow features, and experimental measurement are planes marked) vizualizing 
the wide range of scales in the computational domain and their temporal evolution. 
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measurements in Fig. 26. Figure 24 also provides a visualization of the scale variation, i.e. the very small scales in the 
approaching boundary layer to the large range of scales, at and downstream of, the shear layer reattachment region. 
Contour images of the turbulent structure (visualized via iso-surfaces of vorticity magnitude), shown in Fig. 25, 
provides an alternative three-dimensional visualization of the resolution of fine-scale structure in the approach 
boundary layer, followed by the development of much larger structures, as the flow separates from the corner of the 
backward facing step, to form the free shear layer and recirculation zone. The pronounced effect that the 
recompression shock after re-attachment has on the turbulence structure is also evident in this image. 



Figure 25. Turbulence structure visualization via an instantaneous snapshot 
of iso-surfaces of vorticity magnitude. 

Comparisons of time-averaged flow properties with measurements are given in Fig. 26, where velocity and 
pressure results are compared to available experimental data. Note that an ordinate variables having the superscript 
indicates a coordinate tangential to the ramp surface (e.g. x*) or a coordinate normal to the ramp surface (e.g. y*), 
with all spatial dimensions being normalized by the cavity depth, D. When comparing the shear layer velocity 
profiles with the measured data, it became apparent that the simulations were predicting a slight expansion of the 
shear layer into the cavity. As mentioned previously, the experimentalists translated the cavity ramp during testing to 
prevent this expansion (or possibly compression) in the tests since it was desired to have the shear layer detach under 
constant pressure conditions. In order to concentrate our analysis on the predicted spreading rate of the shear layer, the 
computed results have been shifted by a constant angle to line up the bulk location of the shear layer with the measured 
data. This required a 1.8 degree rotation (relative to the step corner) for the results that utilized the Ducros sensor and a 
2.6 degree rotation for the results that utilized the Larsson sensor. This adjustment of the shear layer alignment may be 
due in part to details of the experimental configuration that were neglected in this effort. The assumption of turbulence 
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homogeneity in the spanwise direction and/or the neglection of any effects that the aerodynamic fences and/or cavity 
sidewalls may have had on the flowfield are possible examples. Once shifted, the LES/RANS shear layer velocity 
profiles agree quite well with the measured spreading rate. The recovery rate of the the boundary layer after shear layer 
reattachment is also captured quite accurately. The mismatch seen in the velocity and pressure profiles normal to the 
inclined surface near the point of reattachment (x*/D=2.7) is due to the slight expansion of the shear layer into the 
cavity as noted above. Predictions of the normalized static pressure along the inclined wall also indicate that the 
simulations predict the shear layer to dip too far into the cavity, as evidenced by the under-prediction of pressure 
upstream of the reattachment and the premature rise in pressure as the shear layer reattaches along the inclined surface. 
Coarse and fine grid solution comparisons were reported in Ref. 43 and the solution was found to be very sensitive to 
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Figure 26. Comparisons of normalized velocity and pressure for FV-SRF with the Ducros (blue long-dash line) 
and Larsson (green short-dash line) with experimental measurements (open triangles). 

grid resolution with the coarse grid solutions yielding shear lay spreading rates that compared poorly with the 
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experimental data. Analysis of the sub-grid modeled turbulent kinetic energy revealed that it was unacceptably large 
relative to the resolved turbulent kinetic energy. 

Further comparisons of the experimental data with the SBP-ASF/PPM-HLLC and SBP-ASFA/VENO-HLLC 
hybrid advection schemes on the 35xl0 6 cell grid are planned but have yet to be completed. Also planned are an 
investigation of modeling assumptions used to reduce the computational cost of the simulations. These include 
geometric effects such as the finite width of the cavity and presence of side fences, grid effects such as the use of non- 
C(0) patches in the spanwise direction at the lip of the backward facing step, and modeling effects such as the use of 
static versus dynamic SGS closure models. 

V. Summary and Conclusions 

This paper provides a study of the components used to form hybrid-advection schemes intended for use in 
Large Eddy and hybrid Reynolds Averaged Navier Stokes/Large Eddy Simulations. The components, consisting of 
non-dissipative advection schemes, dissipative discontinuity-capturing advection schemes, and discontinuity sensors 
used to form the hybrid-advection scheme were implemented in the NASA high-speed propulsion code VULCAN. A 
series of benchmark problems were simulated to test the components against analytic solutions, desirable higher-order 
conservation properties, and experimental data. 

The non-dissipative advection schemes tested were a higher-order finite volume symmetric reconstruction 
flux (FV-SRF) scheme and a higher order summation by parts skew-symmetric ex — split flux (SBP-ASF) scheme. 
The dissipative discontinuity-capturing schemes tested were finite volume, primitive variable reconstruction, single 
quadrature point per face implementations of the k= / 3 MUSCL, PPM, and 4 th -order WENO reconstruction 
schemes. The MUSCL And PPM schemes were tested using the HLLC approximate Riemann solver, and the 4 th -order 
WENO scheme was tested using HLLC and the LLF approximate Riemann solvers. Two benchmark problems were 
computed with all schemes: 1) a grid resolution study of the advection of a 2-D inviscid isentropic vortex, 2) a single 
grid study of the temporal evolution of the inviscid 3-D Taylor-Green vortex. A third benchmark problem was 
computed only using the non-dissipative advection schemes. This problem was the LES of an incompressible, fully- 
developed turbulent flow in a channel. 

The 2-D vortex advection benchmark computations revealed: 

• The MUSCL scheme did a poor job of preserving the vortex on the all three grids. 

• The PPM and 4 th -order WENO schemes did a poor job on the coarse grid and an adequate job on the medium 
and fine grids. 

• Both non-dissipative advection schemes did an adequate job on the coarse and an excellant job on the 
medium and fine grids. 

• Both non-dissipative schemes were shown to be higher-order accurate with the FV-SRF scheme achieving 
near 3 rd -order accuracy and the SBP-ASF scheme achieving design order, i.e. 4 th -order, accuracy. 

• The superior accuracy of the SBP-ASF scheme was due to a much more uniform distribution of error outside 
of the vortex after one pass of the vortex through the computational domain. 

The 3-D Taylor-Green vortex benchmark computations revealed: 

• The dissipative discontinuity-capturing schemes had poor conservation of the mean kinetic energy and mean 
specific thermodynamic entropy (with the 4 th -order WENO scheme using the HLLC approximate Riemann 
solver having the best conservation and MUSCL using the HLLC approximate Riemann solver having the 
worst) emphasizing the need for a non-dissipative scheme when the flow is smooth. 

• The FV-SRF scheme initially dissipated kinetic energy followed by a rapid increase of kinetic energy. 

• The SBP-ASF scheme had superior temporal conservation of the mean kinetic energy while producing a more 
uniform spatial distribution of local kinetic energy compared to the FV-SRF scheme. 

• While both non-dissipative advection schemes schemes have similar mean thermodynamic entropy 
conservation temporal behavior they have very different spatial distributions of local thermodynamic entropy. 
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The grid resolution study of the time averaged behavior of the LES of fully developed turbulent flow in a channel 
revealed: 


• The FV-SRF and SBP-ASF advection schemes provide very similar FES simulation results, with the FV-SRF 
results being marginally better than the SBP-ASF results when compared to DNS data. 

• The level of agreement with the DNS data is quite good for the Reynolds shear stress and the rms streamwise 
velocity fluctuation profiles. 

• The current fine-grid FES results consistently underpredict the DNS values, with a worst case difference of 
less than 5 percent. 

• The current medium-grid LES results show a noticeable improvement over a previous LES effort performed 
on a grid having a similar resolution. 

Sod's 1-D shock tube benchmark problem was computed to examine the ability of three dissipative 
discontinuity-capturing advection schemes to simultaneously capture a rarefaction wave, a contact discontinuity, and a 
shock. A finite volume, primitive-variable reconstruction, single quadrature point per face, implementation of the 
ac = 1/3 MUSCL, PPM, and 4 th -order WENO reconstruction schemes were examined. The MUSCL and PPM 
schemes were tested using the HLLC approximate Riemann solver and the 4 th -order WENO scheme was tested using 
the HLLC and the LLF approximate Riemann solvers to examine the dissipation of the approximate Riemann solver 
on the discontinuity-capturing behavior of the WENO scheme. 

• When using the PPM and the 4 th -order WENO schemes, the most monotone behavior in the vicinity of the 
contact discontinuity was achieved using density, velocity and static temperature as the reconstruction 
variables while density, velocity and static pressure worked equally well for the MUSCL scheme. 

• The LLF approximate Riemann solver was more dissipative than the HLLC scheme with the additional 
dissipation being most noticeable near the leading and trailing edges of the rarefaction wave. 

• The PPM scheme and the 4 th -order WENO scheme gave similar results with the PPM scheme being slightly 
less dissipative than the 4 th -order WENO scheme near the leading and trailing edges of the rarefaction wave. 

Supersonic viscous flow over a cylinder in crossflow was computed to provide a complex temporally 
evolving flow that tests the ability of the hybrid advection scheme to robustly and reliably capture and resolve shock- 
shock and shock-vortex interactions. The discontinuity sensor functions that are key to the success of the hybrid- 
advection scheme were examined to determine their spatial behavior. The Ducros discontinuity sensor function was 
found to force the exclusive use of the non-dissipative scheme, when the sensor had a value of 0, in a small portion of 
the computational domain limited to, the flow upstream of the bow shock, the flow in the cylinder boundary layer, and 
the flow in the cylinder wake. The Larsson discontinuity sensor was found to cause utilization of the non-dissipative 
advection scheme nearly everywhere in the computational domain except in a region encompassing the shocks. The 
improved behavior of the Larsson discontinuity sensor was found to be due to the threshold function and a constraint 
that limited the use of the dissipative discontinuity-capturing advection scheme to regions where the flow was 
undergoing compression. Examination of the instantaneous flow solution discontinuity sensor contours also revealed 
that the overall behavior of the Larsson discontinuity sensor switch was generally good but that sensitivity to the 
threshold function gave rise to fragmentation of the discontinuity sensor contours that could be detrimental to solution 
robustness. The Favre averaged solution computed using the VULCAN code was compared with results obtained from 
a formally 4 th -order accurate finite-difference code (HOFD) on the same grid as well as with HOFD results computed 
on coarser and finer grids. All three hybrid-advection scheme configurations tested in the VULCAN code were found 
to compare well with the time averaged higher-order finite-difference results obtained using the HOFD code. 

A hybrid Reynolds Averaged Navier-Stokes/Large eddy simulation of supersonic flow over a cavity ramp 
configuration typical of those found in high-speed propulsion flow paths was computed using the FV-SRF/PPM-HLLC 
hybrid advection scheme on grids having 8.8xl0 6 and 35xl0 6 cells. The flow was found to be very sensitive to grid 
resolution with the coarse grid solutions having a large fraction of the total turbulent kinetic energy contained in the 


38 

American Institute of Aeronautics and Astronautics 



modeled sub-grid scales. Comparison of the fine grid solutions with experimental data showed that the recovery rate 
of the boundary layer after shear layer re-attachment compared very favorably, and that the predicted free shear layer 
spreading rate compared well with the experiment. There was, however, a disagreement between the simulation and 
experiment with respect to the detachment angle of the free shear layer from the backward facing step lip by as much 
as two and a half degrees. The reason for this difference is currently under investigation. Tests of the effect of the 
discontinuity sensor on the fine grid revealed that the overall results were relatively insensitive to the choice of sensor. 
Further computations are planned to investigate modeling assumptions used to reduce the computational cost of the 
simulations. These effects include geometric effects such as the finite width of the cavity and presence of side fences, 
grid effects such as the use of non-C(O) patches in the spanwise direction at the lip of the backward facing step, and 
modeling effects such as the use of static versus dynamic SGS closure models. Further comparisons of the 
experimental data with the SBP-ASF/PPM-HLLC and SBP-ASF/WENO-F1LLC hybrid advection schemes are also 
planned. 

Overall, although the S BP- ASF low-dissipation advection scheme performed better than the FV-SRF low- 
dissipation advection scheme on the fundamental benchmark tests, the superiority of the SBP-ASF low-dissipation 
advection scheme has yet to be demonstrated on an engineering scale computation and will be the focus of future 
work. With regard to the discontinuity capturing schemes, the PPM scheme and the 4 th -order WENO scheme were 
found to produce very similar results for the cases tested when the same approximate Riemann solver was used. The 
two discontinuity sensors tested were found to produce very different behavior of the hybrid advection blending 
coefficient away from the regions of the flow dominated by vorticity, but to produce similar results in terms of a 
hybrid RANS/LES simulation of an engineering scale flow such as the Settles case. The Larsson discontinuity sensor 
was therefore preferred due to the reduction in computational cost that could be realized by utilizing the expensive 
discontinuity capturing advection scheme only in a a small fraction of the computational domain. The grid sensitivities 
observed in the Settles case warrants further examination of advection schemes that are more accurate than the 4 th - 
order schemes tested. A 6 th -order SBP-ASF scheme / 6 th -order accurate WENO hybrid advection scheme may enable 
such an approach. Further work is also planned to examine additional fundamental benchmark cases such as 
incompressible and compressible isotropic decay of homogeneous turbulence to establish the ability of the candidate 
low-dissipation advection schemes to simulate the turbulent cascade process for the incompressible case, and to test 
the ability of the hybrid advection schemes to robustly capture the eddy shocklets that can spontaneously appear in the 
computational domain for the compressible case. 
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